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ABSTRACT 

A method to rapidly estimate the Fourier power spectrum of a point distribution is presented. 
This method relies on a Taylor expansion of the trigonometric functions. It yields the Fourier 
modes from a number of FFTs, which is controlled by the order N of the expansion and by the 
dimension D of the system. In three dimensions, for the practical value iV = 3, the number 
of FFTs required is 20. 

We apply the method to the measurement of the power spectrum of a periodic point 
distribution that is a local Poisson realization of an underlying stationary field. We derive 
explicit analytic expression for the spectrum, which allows us to quantify — and correct for — 
the biases induced by discreteness and by the truncation of the Taylor expansion, and to bound 
the unknown effects of aliasing of the power spectrum. We show that these aliasing effects 
decrease rapidly with the order N . For = 3, they are expected to be respectively smaller 
than ~ 10"^ and 0.02 at half the Nyquist frequency and at the Nyquist frequency of the grid 
used to perform the FFTs. The only remaining significant source of errors is reduced to the 
unavoidable cosmic/sample variance due to the finite size of the sample. 

The analytical calculations are successfully checked against a cosmological A^-body ex- 
periment. We also consider the initial conditions of this simulation, which correspond to a 
perturbed grid. This allows us to test a case where the local Poisson assumption is incorrect. 
Even in that extreme situation, the third-order Fourier- Taylor estimator behaves well, with 
aliasing effects restrained to at most the percent level at half the Nyquist frequency. 

We also show how to reach arbitrarily large dynamic range in Fourier space {i.e., high 
wavenumber), while keeping statistical errors in control, by appropriately "folding" the parti- 
cle distribution. 

Key words: methods: analytical, data analysis, numerical, statistical, A^-body simulations - 
cosmology: large-scale structure of Universe 



1 INTRODUCTION 

The power spectrum, P{k), represents the primary tool to char- 
acterize the clustering properties of the large scale structure of the 
universe. Most of major constraints on cosmological models and on 
cosmological parameters have been derived from measuring P{k) 
or its Fourier transform, the two-point correlation function. For in- 
stance, the tight constrains derived from WMAP experiment rely on 
measurements of the power spectrum in spherical harmonic space 
(e.g., Dunkley et al., 2008); the most significant results from weak 
tensing analysis come from measurements of the two-point corre- 
lation function of the cosmic shear (e.g., Benjamin et al., 2007; Fu 
et al., 2008); the analysis of the power spectrum of absorption lines 
of lyman-Q forest allowed one to infer drastic constraints on the 
clustering properties of the matter distribution at small scales (e.g.. 
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Croft et al., 1999); and, last but not least, the two-point correla- 
tion function and the power spectrum have been used extensively 
to analyse directly the clustering properties of 2 and 3 dimensional 
galaxy catalogs (e.g., Peebles, 1980; Baumgart & Fry, 1991; Mar- 
tinez, 2008, for a recent general review on the subject). 

To be able to derive predictions from models of large scale 
structure formation, there has been successful attempts to find uni- 
versal dynamical laws, partly phenomenological, that lead to semi- 
analytical expressions of the non linear power spectrum (or the two- 
point correlation function) of the matter distribution. Among them, 
one can cite the nonlinear ansatz of Hamilton et al. (1991), later 
improved by Peacock & Dodds (1996, see also Smith et al., 2003). 
Such a non-linear ansatz has been used to constrain models against 
observations, particularly in weak tensing surveys (e.g., Benjamin 
et al. 2007; Fu et al. 2008). Another well known phenomenological 
description is the so called halo model, which proposes not only 
some insights on the clustering properties of the dark matter dis- 
tribution, but also of the galaxy distribution itself (see, e.g.. Ma & 
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Fry, 2000; Peacock & Smith, 2000; Seljak, 2000; Scoccimarro et al. 
2001; see Cooray & Sheth 2002 for an extensive review). Again, the 
measurement of two-point statistics in the galaxy distribution was 
used to infer important constraints on the halo model parameters 
(e.g., Abazajian et al. 2005). 

With the advent of very large modern surveys, "precision cos- 
mology" has become a reality: the accuracy of the observations 
have caught up to the accuracy of the predictions. These need to be 
more and more precise to constrain cunent models of large scale 
structure formation, for instance the fiducial ACDM (Cold Dark 
Matter) model. It is therefore now crucial to obtain very fine esti- 
mates of statistics in simulations with good control of the system- 
atic errors in order to be able to fine tune non linear ansatz such that 
of Hamilton et al. (1991) or details on the set up of the halo model. 

In this work, we concentrate on the problem of measuring as 
accurately as possible the Fourier modes of a distribution of points 
such as those coming from an A'^ body simulation, with a partic- 
ular emphasis on the power spectrum. The traditional method for 
measuring the Fourier modes, 5k, consists in assigning the point 
distribution to a periodic grid with some interpolation method and 
then computing 5k with a Fast Fourier Transform (FFT) technique. 
However, the introduction of a grid, combined with the correspond- 
ing interpolation, induces two effects: damping of the modes at 
large k due to the convolution involved in the interpolation, and 
effects of aliasing due to the finite resolution of the grid (e.g., Hock- 
ney & Eastwood, 1988). In addition, the discrete nature of the par- 
ticle distribution induces some systematic effects, but these latter 
can be straightforwardly accounted for if the distribution of parti- 
cles is a local Poisson process (e.g., Peebles, 1980). While the bias 
induced by the interpolation method can also be easily corrected 
for, the effects of aliasing are more difficult to control. This has 
been for instance illustrated well by Jing (2005), who proposed to 
correct for aliasing using an iterative method, based on an ansatz 
that assumes that the power spectrum behaves like a power law 
at large k. However this method, although efficient for cosmologi- 
cal power spectra which behave close to power laws, is not free of 
biases in general. An alternative route involves using appropriate 
interpolation functions, which are by construction meant to reduce 
the effects of aliasing as much as possible. This is for example the 
case of the Daubechies wavelets (Daubechies, 1988), as proposed 
by Cui et al. (2008). While these interpolating functions are power- 
ful, there are still some significant residuals at the few percent level, 
and Cui et al. (2008) do not provide a rigorous way to quantify and 
correct them. 

The aim of this paper is more ambitious than Jing (2005) and 
Cui et al. (2008): we want to be able to measure the power spectrum 
from a simulation at an arbitrary level of accuracy, with rigorous 
control of the biases and the residuals due to aliasing. Of course, 
the higher the required level of accuracy, the larger the computa- 
tional cost. Furthermore, even though we shall be able to measure 
the power spectrum extremely accurately from a given simulation, 
it does not mean that the power spectrum of the underlying cosmol- 
ogy will be estimated fairly: a statistical error — cosmic variance 
— arising from the finite number of available modes given the finite 
size of the simulated volume will still be present (e.g., Feldman, 
Kaiser & Peacock, 1994; Scoccimarro, Zaldarriaga & Hui, 1999; 
Szapudi, 2001; Bernardeau et al. 2002 for a review). 

The method we propose is inspired by Anderson & Dahleh 
(1996). It is based on the fact that the Taylor series expansion of 
trigonomic functions, sin(a;) and cos(a;), converges very rapidly. 
This will allow us to compute Fourier modes efficiently with a num- 
ber of FFT depending on the order A'^ of the Taylor series expan- 



sion used. That number will control the effects of aliasing as well 
as the biases on the rough estimator, which can be corrected for in 
the case of a local Poisson realization of a stationary random field. 
We shall write explicit analytical expressions for the power spec- 
trum and propose an unbiased estimator that will be tested in detail 
against a controlled A^-body experiment. 

This paper is organized as follows. First we describe what we 
call the Fourier-Taylor transform and its practical implementation 
(§|2j. Then, we construct a rough estimator of the power spectrum 
from it, and perform analytical calculation of its ensemble average 
by assuming local Poisson sampling of a stationary random field 
(§[3]l. This section is supplemented with Appendix A, which dis- 
cusses some subtle differences between the unconstrained versus 
the constrained ensemble average, and Appendix B, which details 
some useful analytic expressions of various quantities occurring in 
the calculations. We study the biases on the rough estimator of the 
power spectrum, which can be easily corrected for, as well as the 
unknown residuals due to aliasing, which are controlled by the or- 
der of the Taylor expansion. The analytic results are then validated 
in a CDM GADGET A^-body simulation (§|4]l. We study two config- 
urations, the final stage of the simulation, which should agree very 
well with the assumption of local Poisson sampling and stationar- 
ity, and the initial conditions, corresponding to a slightly perturbed 
grid. This section is supplemented with Appendix C, which details 
the calculation of the power spectrum of a randomly perturbed grid. 
In § [5] we show how to cover all the available dynamic range in 
Fourier space while keeping control of the errors, by appropriate 
foldings of the particle distribution. Finally, section[6]briefly sum- 
marizes the results of this paper 



2 THE FOURIER-TAYLOR TRANSFORM 

We begin with a discrete distribution of points Xi with weights 
(masses) Mi, i = 1, • • • , A^p, where each Xi is potentially a D- 
dimensional vector. The equivalent perturbed density field is 

1 

p{x) = —^Ah5D{x-x,) (1) 

i = l 

where 5d{x) is the D-dimensional Dirac delta function. The 
Fourier transform of this distribution is 

5k = J d^x p{x) e' '^-^ 

= -^^M,exp(Jfc-2;0> (2) 

^ i 

where fc is a D-dimensional wavevector and the imaginary unit is 
= — 1 (we use lower case i as an integer index). Since the num- 
ber of dimensions is assumed to be arbitrary the operator "■" is the 
scalar product. The direct calculation of the sum in equation Q 
is a very slow, an A^^; x A^p process, where A^^ is the number of 
sought wavenumbers. 

To speed up the calculation, one can choose a homogeneous 
cubic grid of a certain size covering the volume occupied by all the 
points, Ng (in 3D, Ng x Ng x 7Vg)Qand define the function J(i) 
giving the (vector) integer position of the cell containing object i. 
Then equation l|2j becomes 



^ The following calculations can be easily generalized to a rectangular grid. 
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^^^i^E Af.exp{7fc.[i + A(i)]}, 

with 

A(i) = x,- J{i). 



(3) 



(4) 



The quantity A(i) (or each coordinate of it in more than ID) is 
bounded in [—1/2, l/2[. To simplify the expressions, it is assumed 
without any loss of generality that the size of a cell of the grid is 
unity, Ag = 1. The Fourier- Taylor expansion of order A'^ is then 



r(iV) 



iV 

j n — i\J{i)=j 

Such an expansion is expected to converge very quickly as a conse- 
quence of properties of trigonometric functions, sin(2;) and cos(a::). 

With D the number of dimensions, and the vector A = 
(Ai, ■ • • , Ad), from the multinomial theorem, 

(fc-A)" = 



V 

91 H I-ijd=" 

x(fciAi^''i 



X qo'. 
X ••• X (/cdAd)"". 
Equation (O can thus be rewritten 

N 

1 



(6) 



X go! 



71—0 giH Vqu—"!^ 

xkl^ X ■■■ X k']^F{k,q) 

with 

F{k, q) = FTK] = Z'^^^) ^^pI-^'^ ■ ^1' 
where FT is the Fourier operator and 
M,(j)= > ; x[Ai(i)]^^ x...x[AdW]'". 



qi\ X 



E 

i\j(i)=] 



(7) 
(8) 

(9) 



This defines the Fourier- Taylor algorithm: the approximate direct 
Fourier transform is reduced to 

(i) the calculation of the moments on the the real space 
grid, equation (|9}; 

(ii) their Fourier transform, equation ([8](, which can be per- 
formed with usual FFT algorithms; 

(iii) their summation with the appropriate weights, equation (O. 

Note importantly that periodicity was not assumed in this calcu- 
lation, and that the values of (each coordinate of) k available are 
theoretically any multiples of 2n/Ng, as a simple consequence of 
the periodicity of the function F{k, q) in k space. However, the ac- 
curacy of the Taylor expansion is controlled by the magnitude of 
k ■ A(i), and thus worsens with larger k. As a result, we shall re- 
strict at present time to the naturally available range of values of 
(each coordinate of) k, [— fcny,fcny], where fcny = tt corresponds 
to the Nyquist frequency defined by the grid, which is a priori un- 
related to the distribution of points. We shall explain in §[5]how to 
extend the algorithm to have access to arbitrary values of k, while 
maintaining the errors on the Taylor expansion bounded. 

It is crucial to point out a few features of this calculation. First, 
this is very specifically the Fourier transform of a point distribution, 
which is not precisely equivalent to the transform of an irregularly- 
sampled continuous function (which requires the further specifica- 
tion of an interpolation scheme). Second, if we restrict the calcu- 
lation to a finite set of wavenumbers k, there is no unique inverse. 



Table 1. Number of Fourier transforms, A'^ppT and the rough estimate of 
the enor, E(D, N). in the worse case (at the Nyquist frequency) according 
to equation jilt , as functions of the considered order A'^ of the Fourier- 
Taylor expansion and the number of dimensions D of the system. 
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and we cannot recover the real-space distribution from the Fourier 
transform. Finally, the lowest-order (TV = 0) version of this cal- 
culation is equivalent to nearest-grid-point interpolation to the Nj^ 
grid. 

The algorithm now scales in three dimensions like C[A''fft x 
Nl log TVg] + OfiVpFT X iVp], for accessing Nu ~ iV"! wavenum- 



bers, where A^fpt represents the number of Fourier transforms in- 
volved in the calculation. If one assumes A*'p A'^, this method is 
much faster than the direct summation approach if A^ppx <C N^- 
The parameter A'ppt is given by 



iVpFT = ^ ^ 



(n + L>)! 
D\ n\ 



(10) 

1=0 91 H h'jr>=i 

Table [T] gives the corresponding numbers for D = 1, 2 and 3. The 
accuracy of the approximation is dictated by the magnitude of the 
next order correction in equation (|5]l, [fc- A(i)]^"'"i / (A^+1) !. Errors 
become more significant at the Nyquist frequency of the grid and 
for A(i) ~ 1/2. At first sight, control of the error is thus given by 
the condition 



E[D, N) = (7rLi/2)™ + 7(iV + 1) 



)! <:, 



(11) 



where e is a small parameter, but of course the actual error depends 
on the spectral properties of the system considered. While order 
A'' = 10 is enough to have e < 10"'' for D = 1, we need A'^ = 15 
and A'^ = 20 for D = 2 and D = 3, respectively. For this level of 
accuracy, the computational cost becomes increasingly prohibitive 
for increasing value of D due to the large number of Fourier trans- 
forms required to perform the calculations. This makes the Fourier- 
Taylor approximation mainly attractive for D = 1 if one aims to 
estimate 5^ accurately for any value of k. However, the goal here is 
not to have an accurate measurement of 5k but rather of its power 
spectrum. Let us now investigate how the Fourier-Taylor method 
behaves for this latter quantity. 



3 THE FOURIER-TAYLOR POWER SPECTRUM 



This section is divided into two parts. In §[3TT] we compute the en- 
semble average of the naively-defined rough Fourier-Taylor power 
spectrum, assuming that the point process under consideration is a 
local Poisson realization of a stationary random field, periodic over 
the grid used to run the Fourier- Taylor algorithm. For instance we 
shall recover a well-known result for nearest grid point interpola- 
tion (NGP), which corresponds to the zeroth-order Taylor expan- 
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sion. In § 13.21 we analyse the various biases in the rough Fourier- 
Taylor estimator, namely the shot noise of the particles which can 
be subtracted off, the bias due to the interpolation method (e.g., the 
famous sinc^ biasing from NGP interpolation), which can be easily 
corrected for, and the residual due to aliasing. The calculations here 
do not necessarily assume isotropy of the underlying random pro- 
cess, but the 3D analyses use angular averages, which make sense 
only if isotropy applies. 



3.1 Ensemble average for a stationary point process 

In what follows, we assume that the catalog is a set of particles of 
equal weights. Mi = 1. A naive estimate of the power spectrum of 
order N can be written 



(12) 



where 



K S^'S'-V = > ' exp[/fc ■ U - j')] V X 



j2 x(N)j^{N) _ 

E E 



n.m — 



X 



[ik-A{i)r[-ik-A{i')r, (13) 



iVp — (Np), and (• • ■) stands for the ensemble average over many 
realizations. Note thus that in the following calculation, we shall 
allow the number of objects A^p in the catalog to fluctuate. Also, 
isotropy is not yet assumed here: k is still a vector in equation il2l . 
Setting 



Vn{k) 



and 



1/2 



(14) 



1/2 



1/2 



Vn,m{k,j-j') = / C(j-j'+A- A') X 

J -1/2 

x(A:-A)" (fc- A')'"d°Ad°A', (15) 

where ^(a;) is the two-point correlation function assumed to be in- 
variant by translation, ensemble averaging equation l ll3t reads 



y i 

ri,m— 



{iVp Vn,+m{k) + 



+ N'Y, exp[/fc ■ (i " /)] J'n,™(fc, j - j')}, (16) 

3,3' 

where S-olk) is the Dirac delta function and N is the average num- 
ber of particles per cell, 

N = iVp/iV|. (17) 

This calculation can be derived quite easily following the micro- 
cells formalism of Peebles (1980), as explained in Appendix lAl I 
Let us assume periodicity over the grid and decompose the func- 
tion ^(r) in Fourier modes. 



C(0 



P(0exp(-7rr), 



(18) 



where P{1) is in fact the sought power spectrum. Notice that the 
sum {Tsj is infinite because the system is not necessarily band- 



limited. Then, 

exp[7fc . (j - j')] Un,m{k,j - j') = 

j.j' 

Y Pil) eMHk-l)-U-j')]^n,mil,k), (19) 

with 



{l,k) 



/ exp[-7 / ■ (A - A')] X 
J-1/2 

x(/fc- A)" {-Ik-A')"'dA dA' 



(20) 



The sums over j and j' cancel unless I = k + 2ttM, where M is 
an arbitrary (vector) integer. Thus 

jn-m ^ exp[/A; ■ (j - j')]iyn.m{k,j - f) = 
j.j' 

Ng Y + 27rM)K„,,„(fc + 27rM, k). (21) 
Notice that 

Kn,m{k + 2tv M,k) = K„(fc,Af) X K.m{k,M), (22) 

with 

K^{k,M)= / exp[-J(fc + 27rM).A] (/fc- A)"d°A.(23) 

J-l/2 

Details of the calculation of the (real) number K„(fc, M) are given 
in Appendix B. We thus obtain the simple expression 



P'"'{k) = 5u{k) + ^WN{k) + 

JVp 

+ ^P(fc + 27rM) Tl{k,M), 



with 



rN{k,M) = J2 



Hn{k,M) 



n = 
JV 



WN{k) = y -^—:iy,n+nik). 



(24) 



(25) 



(26) 



It is easy to check from equations l ll4b and l l23t that 

lim Tjv(fc,M) = (5d(M), lim W^jv(fc) = 1, (27) 

N—*oo iV— »oo 

as expected: the Fourier- Taylor approximation tends to the exact 
solution, P{k) + 1/ Np, when the order N oo. 

Using the analytical calculations presented in Appendix B, we 
recover for A'^ = (corresponding to NGP interpolation), the well- 
known result {e.g., Jing, 2005) 

P^''\k) = 5D(fc) + i +^P(A; + 27rM) X 

^ u 

X Y\ [sin(fc,/2)/(fc,/2 + 7rM,)]'. (28) 

9=1, ■■•,D 

The bias on the power spectrum introduced by the NGP interpo- 
lation corresponds, in the ensemble average sense, to the Fourier 
transform of the square top-hat function modulo aliasing of the 
power spectrum in Fourier space, and can be corrected for straight- 
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Figure 1. Tjv(fc, M) as a function of k/k^y + 2M in the ID case, for 
various values of the order TV as indicated on the panel. The calculation 
has been performed using A^g = 128, but the results would not change 
significantly for other values of Ng . 



forwardly in the band-limited case (where only M = con- 
tributes)!; 



3.2 Analysis of biases and residuals due to aliasing 



The problem with equation l l24t is the sum over M, as it corre- 
sponds to foldings of Fourier modes at values of k we do not have 
access to for a given grid size. The higher the order considered, the 
smaller the effect of this aliasing, by construction of the Fourier- 
Taylor method. This is illustrated in ID by Fig. [T] which shows 
the function Tjv(fc, A/) as a function of k/kny + 2Af, for vari- 
ous orders of the Taylor expansion. One can see that convergence 
towards the exact solution (equation I27t is rather fast. The effect 
of aliasing is largest when approaching Nyquist frequency, a nat- 
ural property of the Fourier-Taylor method which is an expansion 
around k = 00 While it is not possible to correct for aliasing of 
the power spectrum without additional strong prior assumptions, it 
is possible to estimate a bound on the systematic error it induces. 
For instance, let us assume that outside the range [ /Cny^^ny] (for 
each coordinate of the wavenumber in more than ID) the power 
spectrum is bounded by a value Pmax: 



P{k) ^ Pmax, outside Nyquist range. 



(29) 



^ Note as well that for pure white noise, P{k) =constant, NGP is unbiased 
since ^^^l/(fe,/2 + 7rAfq)2 = 1/ sin2(fe,/2). 

^ Note that a similar polynomial expansion (or an expansion on an 
appropriate basis of functions) to the Taylor expansion but minimiz- 
ing in a global way the effects of the foldings would be more op- 
timal. One could for instance adapt methods used in the NFFT al- 
gorithm (www-user . tu-chemnit z . de/~pott s/ nfft/doc. php, 
see, e.g.. Potts, Steidl & Tasche 2001), potentially more efficient than the 
Fourier-Taylor transform. The advantage of this latter method presented 
here is the simple analytic control of all the biases. 



This obtains, for example, in the common case of a falling spectrum 
at high wavenumbers. Then, noticing that (whatever the number of 
dimensions) 



J2'^N{k,M) ^WN{k), 

M 



(30) 



we obtain the following bound (omitting the additional trivial addi- 
tional term at fc = 0) 

pW(fc)-iy^(fc)/jVp 

where the positive residual function _Rjv (fc) is given by 
WN{k) ^ 



P{k)+P^^^RN{k), (31) 



Piv(fc) 



(32) 



Equation ( I31t defines a range for the estimation of the unbiased 
power spectrum with the weak assumption given by equation l |29t . 

The residual 7?jv(A;) estimates the influence of the foldings 
of the (unknown) power spectrum at wavenumbers outside the 
Nyquist domain defined by the grid. It is expected to decrease 
rapidly with order A'^, since Wjv(fc) ~ T^(fc,0) when \k\/kny < 
1: at leading order in k/kny, after simple algebraic calculations, 
one finds, for N ^ 1, 



T^(fc,0) 

2(-l)^/2 



N\{N + 2) 
2(-i)(^+i)/2 
(iV+1)! 



(33) 

UN+2{k), N even (34) 

m+i{k), iVodd. (35) 
only for Tg(fc,0), while 



Equation i34\ remains valid for A*' 

Wo{k) = 1. 

To illustrate quantitatively these results in the D — 3 case. 
Figure |2] shows the angular averages of the biases T^(A;,0) — 1, 
M^jv(fc) — 1 (left panel) and the residual i?jv(/c) as functions of |fc| 
for various values of A''. More specifically, and this will be used fur- 
ther for 3D measurements, one estimates for each integer wavevec- 
tor kNg/{2n) the following quantity 



fc(fc) = E I 



( kNg l\ 
\ 2tt 2) 



(36) 



where E{x) is the integer part of x. Then, the angular average of 
quantity A(k) for integer wavenumber modulus k is given by 



C{k) 



J2 



fcjfe(fe)=fc 



where the count C{k) is 

c{k)^ 1- 

k\k(k) = k 



(37) 



(38) 



This gives the number of integer wavenumbers verifying k{k) — k. 
Note that angular averages make sense only if one assumes statis- 
tical isotropy, which then means that the power spectrum depends 
only on \k\. This is theoretically the case of the cosmological ran- 
dom fields considered in this work. 

The way the angular average is performed here is very rough, 
and itself introduces some biases with respect to the estimate of the 
true angular average of the power spectrum. Implementation of a 
better angular averaging procedure would be quite straightforward 
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Figure 2. Biases on the Fourier Taylor spectrum, |T^(fc, 0) — 1|, 
|VKjv(fc) — 1|, and residual function Rp{(k) due to aliasing as functions 
of fc/fcny, after angular average as explained in the main text. The calcu- 
lation assumes a grid with Nf, = 128, but the results would not change 
significantly for other large enough values of Nf, . Each curve corresponds 
to a value of the Taylor expansion order, A'^, ranging from A*^ = to A'^ = 6, 
from top to bottom. On the top panel, the solid, dotted and thick grey curves 
correspond to |T^(A:, 0) — 1|, \ Wis[{k) — 1| and leading order expression 
)35K respectively. The case Af = is omitted, for clarity. In this latter case, 
we have W(i{k) — 1 = and |T|^{fc, 0) — 1| is of the same order (but 
slightly different) of what is obtained at first order, A^ = 1. 



and easy to propagate in the analytic calculations, but would not 
serve the purpose of this paper, so we leave it for future work0 

Figure [2] shows that the residual function i?jv(fc) decreases 
rapidly with . It is supplemented with Table ID which provides 



* See, e.g., ScoccimaiTO et al. (1998), for a better handling of angular av- 
erages. 



Table 2. Numerical estimate in the 3 dimensional case of the residual func- 
tion i?,jv(fe) after angular average as explained in the text. The calculation 
has been performed with Afg = 128 but the results should not change sig- 
nificantly for other values of A^g as long as they are not too small. The first 
column indicates the order A^ of the Taylor expansion, while the second and 
the third one give iJjv(^ny/2) and Rpf{kny), respectively. 



N 


fljv(fcny/2) 


ijjv(fcny) 





0.22 


1.3 


1 


0.011 


0.17 


2 


1.3 10-'^ 


0.055 


3 


5.6 10^5 


0.018 


4 


2.0 lO"*' 


2.2 10-3 


5 


5.2 10"* 


2.6 10"* 


6 


1.1 10"^ 


2.0 10-^ 



. One can already see 



numerical values for k = fcny/2 and k — 
the virtue of the Fourier-Taylor method: going to higher order re- 
duces these otherwise uncontrollable effects, which are clearly not 
negligible at all for the traditional NGP method (A'^ = 0) where 
the residual is of order unity at the Nyquist frequency and still 
about 20 percent at half the Nyquist frequency, whereas the third- 
order Fourier-Taylor correction reduces it to about 2 percent and 
6 X 10^^, respectively. 

To understand better the scaling of the functions T^(fc, 0) — 1 
and Wjv(fc) — 1 with k, one can perform the integral of equa- 
tion il4\ in a sphere instead of a cube. For Z) = 3, it reads 



3 / 3 X'^/a 



(47r) 



7r|fc| 



(39) 



2(iV+l)(Ar + 3) 

This approximation is not accurate enough for practical calcula- 
tions, but allows one to see that Wjv(fc) — 1 and T%{k,0) — 1 
scale as for A'' and + 1, A'^ even. For instance, the second 

and third order scale the same way with as can been seen on top 
panel of Fig. [2] Note interestingly the bending of the bias for odd 
N observed when one approaches the Nyquist frequency: this fol- 
lows from the fact that the sine function cancels at fc = fcny, unlike 
the cosine function. Therefore, if one wants to use the brute-force 
Fourier-Taylor expansion without a bias correction to the power 
spectrum, it is better to perform it at odd orders. 



4 VALIDATION: TESTS ON AN W-BODY SIMULATION 

In this section, we validate the analytic results derived just above 
by performing measurements in a CDM A^-body simulation. The 
simulation is described in § 14.11 where we also explain how we 
estimate statistical errors on the measurement of the power spec- 
trum. In § 14.21 biases on the Fourier-Taylor rough estimator are 
measured and compared to the theoretical predictions in the final 
stage of the simulation, which should agree well with the assump- 
tions of local Poisson sampling of a stationary random field used in 
the previous section. In 14.31 we consider the initial conditions of 
the simulation, which correspond to a slightly perturbed grid and 
therefore strongly deviate from local Poisson behavior. Finally, in 
§ 14.41 we propose a nearly unbiased estimator that should work in 
all the cases, disregarding the aliasing effects on the power spec- 
trum, which are controlled by the order of the Fourier-Taylor ex- 
pansion as well as the ratio of fc/fcny where fcny is the Nyquist 
frequency of the grid used to perform the calculations. 
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Figure 3. A thin slice, L/128 tliick, extracted from tlie initial conditions and the final stage of our 128'' CDM GADGET simulation. On the left panel, 
because of the strong deviations from local Poisson behavior brought by the grid pattern, the validity of our calculations for the power spectrum biases are 
questionable. On the right panel, the effects of the grid are much less present, although still visible in underdense regions: they should have much less impact 
on the measurements. 



4. 1 The CDM G AD GE T sample 



We now validate the above results using a Cold Dark Matter 
(CDM) TV-body simulation performed with the publicly available 
tree code GADGET (Springel, Yoshida & White, 2001), as shown 
in Fig. [3l The parameters of this simulation are the following. It 
uses A^p = 128^ particles in a periodic box of size L = 50 
Mpc, where h — iyo/(100km/s/Mpc) — 0.7. The cosmological 
parameters are matter density = 0.3 and cosmological con- 
stant =0.7. The linear variance erg of the density fluctuations 
in a sphere of radius 8 Mpc extrapolated to present time was 
taken to be erg — 0.92. Finally, the softening length e was chosen 
to be l/20th of the mean interparticle distance. The initial condi- 
tions were generated using the Graphic package of Berstchinger 
(2001), with a transfer function given by Bardeen et al. (1986) (no 
baryons). This package basically allows one to perturb an initially 
homogeneous particle distribution using the Zel'dovich approxima- 
tion (Zel'dovich, 1970). For the initial pattern, we chose to put the 
particles on a regular grid (see left panel of Fig.jSj. Our test sample 
is somewhat small compared to contemporary numerical experi- 
ments, but was chosen such that we could perform "exact" calcula- 
tions of power spectra in a reasonable amount of time. By exact, we 
mean the 20th order Fourier-Taylor expansion using an A'^g = 128 
grid. To be able to probe the highly nonlinear regime well and to 
make sure that the system has properly relaxed to a locally Pois- 
sonian and isotropic stage in collapsed objects (see right panel of 
Fig. [3}, we purposely used a small box size. 

The "exact" power spectrum of the particle distribution in this 
sample is shown in Fig. |4l both for the initial conditions and the 
final stage of the simulation. The smooth curves correspond to the 
nonlinear ansatz of Hamilton et al. (1991) using the formula of Pea- 
cock & Dodds (1996), shown here for reference. The error bars 
represented by the gray shadded areas correspond to the following 



10° 



10" 



Present time 




CDM GADGET sample 



10 

kL/27T 



Figure 4. The power spectrum measured in the initial conditions (lower 
curves) and the last snapshot, corresponding to the present time (upper 
curves), of our 128^ CDM GADGET sample. The symbols correspond to 
the "exact" measurement with the Fourier-Taylor method of order 20. A 
shaded region is superposed on them. This represents uncertainties on the 
measurements as computed from equation {40\ . The smooth curve gives the 
non-linear ansatz of Peacock & Dodds (1996). The horizontal dotted line 
corresponds to the white noise level of the particle distribution. Note that 
for the measurement at the present time, a connection for white noise was 
performed, but not for the initial conditions measurement. The bending of 
the power spectrum at high k in the latter case is due to the Manning filtering 
performed in the Graphic package of Berstchinger (2001). 
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self-consistent estimate of tfie statistical error: 

2 



^-frough (^) 
-frough (^) 



C(k)[C{k) - 1] 



E 



(<5fc5-fc)'"C(fc)P,Lgh(fc) 



where 

-frough {^) 



C(fe) ^ 

k/k(k) = k 



SkS- 



(40) 



(41) 



is the rough power spectrum measured from the distribution of par- 
ticles. From now on we omit the cumbersome tilde on P and bar on 
k as we shall only consider angular averages until the end of this 
section. In the framework of an isotropic, stationary local Poisson 
process discussed in §[3] recall that after ensemble averaging over 
many realizations 



{Prough(fc)) = P{k) + — 



(42) 



where P{k) is the underlying power spectrum. We noticed that 
the error given by equation i40\ is well-approximated by the well- 
known result obtained for a random Gaussian field (e.g., Feldman, 
Kaiser & Peacock, 1994) 



C{k)' 



(43) 



which translates, for the desired shot-noise-corrected power spec- 
trum to 



^Py 



dk) 



1 + 



+ 



NpP{k) N^P^{k) 



(44) 



Here, C{k) represents the number of .statistically-independent 
available wavenumbers, hence the missing factor of two compared 
to the usual case, since symmetries in fc-space are already taken 
into account. We are a bit puzzled by the very good agreement be- 
tween equation J40t and l l43t . as we would expect non-Gaussian 
contributions to the error in equation i40l due to nonlinear cou- 
pling generated by the dynamics. It seems that these couplings are 
rather small, as already noticed by Rimes & Hamilton (2006) and 
Hamilton, Rimes & Scoccimarro (2006). Of course, we know for 
sure that the error given by equation l |40| l and ( 143 1 underestimate 
the true value in general, that would be obtained from the dipersion 
over many realizations of our simulation (Scoccimarro, Zaldarriaga 
& Hui, 1999). However, the realistic calculation of such an error 
requires prior knowledge of the bispectrum and the trispectrum. 
Moreover, eqs. l l40t and l l43t are sufficient to prove the points dis- 
cussed in the analyses of this paper and to provide a rough estimate 
of errors in a simple and self-consistent way, as can be provided by 
the numerical package we propoself] 

Due to the very small size of the box, the agreement between 
the smooth curves and the measurement on Fig.|4]is not very good 
at large scales (small fc), where few individual modes are available; 



^ For further discussion of statistical errors on the power spectrum, in pai- 
ticulai' some possible improvements of equation )40t for a self-consistent 
calculation of the errors and the covariance matrix of the measured power 
spectrum, see Hamilton, Rimes & Scoccimarro (2006). 



this effect is even worse at the final stage of the simulation be- 
cause of the nonlinear coupling at scales close to the simulation box 
size. Note importantly that the choice of a regular pattern combined 
with the Zel'dovich approximation for the initial particle distribu- 
tion has a non-trivial influence on the evolution of individual modes 
(Marcos et al., 2006; Joyce & Marcos, 2007a,b; see also Crocce, 
Pueblas & Scoccimarro, 2006, who discuss transients coming from 
using the Zel'dovich approximation). Note also that the damping 
of the power spectrum at large k measured in the initial conditions 
is not due to any interpolation effect — as we have access here to 
an "exact" measurement — but to the Hanning filtering performed 
in Graphic (see Bertschinger, 2001). Finally, while a shot-noise 
correction was performed on the P{k) obtained in the final stage 
of the simulation, i.e., a term was subtracted from the rough 

measurement, we reiterate that it does not apply to the initial stage. 
In this case, it is more appropriate to perform no correction as we 
are in the situation of a perturbed grid pattern (e.g., Joyce & Mar- 
cos, 2007a)Q 

4.2 The ideal situation: final, relaxed stage 

Figure[5]shows the measured bias on the shot-noise-corrected mea- 
sured Fourier-Taylor power spectrum of order A'^, for various values 
of iV; 



p(^\k)-WNik)/Np^P{k) 
P{k) 



(45) 



If the analytic calculations of ij |3.2l applv. we should have 

Pm 



ll{k,0)<ib{k)<,'y%{k,0) + 



P{k) 



[VKiv(fc)-T^(fc,0)].(46) 



The lower and upper bound are represented by the continuous and 



the dotted line respectively. We assume Pn, 



- P(fcny). As ex- 
pected, within the statistical errors defined by equation ( I40t . the 
measurements represented by the symbols indeed lie in between 
these two curves. Since the true power spectrum is a strongly de- 
creasing function of k in the CDM cosmology, the effects of alias- 
ing of the power spectrum are overestimated by the dotted curve, 
so the symbols are much closer to the solid curve than to the dot- 
ted one. In fact, they overlay quite well on the solid curve when k 
is small enough compared to k^y, as expected. Note that the near- 
est grid point interpolation (upper left panel) is still very signifi- 
cantly contaminated by aliasing, while this effect decreases rapidly 
with the order A'^ as shown in previous section. The bias also be- 
comes smaller and smaller with N, but it is really worth correcting 
for, since the theoretical predictions match rather well the measure- 
ments. Note that at high order, TV ^ 4, the symbols no longer lie be- 
tween the dotted and solid curve, but this is nonetheless well within 
the statistical errors represented as the shaded region. Moreover, the 
system still deviates locally from a pure random stationary pattern 
(right panel of Fig. |3) so one cannot expect the theoretical bound 
given by equation l l46t to remain valid at such a level of accuracy. 

4.3 A poor situation: a perturbed grid as in the initial stage 

While the final stage of our GADGET sample is well within the 
framework of the assumptions of local Poisson sampling of a sta- 
tionary random process, this is not the case for the initial con- 
ditions, which correspond to a slightly perturbed grid pattern, as 

^ This can be easily checked by analysing equation 1501 which gives the 
power spectrum of a perturbed grid using the Zel'dovich approximation. 



Accurate estimators of power spectra in N-body simulations 9 




0.4 0.6 

k/k„. 



0.3 



0.1 



0.0 
0.0 



order 3 

1 


•1 

■ o/ '. 

■ o/ 
■o/ 


i 


.'o/ J 
•0/ 

o/ 

'7 



0.4 0.6 



0.0030 
0.0025 




0.0000 




ooo Measurement 

Error 
Theory 

Theory + P(k„,) folding 



CDM GADGET sample, present time 



Figure 5. The measured bias, b(k) (eguation l46t . on the rough estimator of the power spectrum, p(^) (fc), as a function of k for our CDM GADGET sample 
(symbols). Each panel corresponds to a value of the order A'^, as indicated. The solid and dotted curves represent a theoretical lower and upper bound, 
respectively, in between which the symbols should lie requation l46l with Pmax = f (fcny)], in the framework of §[5] — local Poisson sampling of an isotropic 
stationary random field. The shaded region represents the statistical eiTor computed from equation (40). 



10 Colombi, Jqffe, Novikov & Pichon 




0.10 
0.08 
0.06 
0.04 
0.02 
0.00 
-0.02 



-0.04 



order 4 


• O 




; O " 




; o " 
■' o 




' 







0.0 0.2 



0.4 0.6 



0.8 1.0 




ooo Measurement 
Error 

Theory 

Theory + while noise folding 



CDM GADGET sample, initial conditions 



Figure 6. Same as Fig. [5] but the function displayed is now g{k) (eQuation l47t as measured in the initial conditions of our GADGET sample, which correspond to 
a perturbed grid. The soUd lines on each panel still show function 7^ (fc, 0), while the dotted line gives 7^ (A:, 0) -I- (f rough {k)N'p)~^[Wj^{k) — (k, 0)] 
as discussed in the text. 
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shown on left panel of Fig. [3] Fig.[6]is the same as Fig.|5] but now 
the function measured is 



P'^Hfc) - Plough (fc) 

-frough (^) 



(47) 



In other words, no correction for the shot-noise of the particles is 
performed since it does not make sense in that case. 

Before pushing this analysis further, we have to understand 
what is the particle sample at hand. The Graphic package per- 
turbs an initial grid pattern with a Gaussian random displacement 
field, Viq), supposed to be stationary, isotropic and curl-free. We 
assume, to simplify the discussion which follows, that the grid of 
particles is the same as the grid used to perform fast Fourier trans- 
forms, A'^p — Ng. This is the case for Fig.[6] Remember finally that 
lengths are expressed in units of the size of the a cell of the grid. 

The perturbed position of the particle is given by 



X = s + q + P{q), 



(48) 



where q is an integer vector for which each coordinate is in the 
range [0, A^g — 1], and s is a small constant offset, for which each 
coordinate is in [0, 1[. Usually, s = or s = (0.5, 0.5, 0.5). The 
Fourier mode of the perturbed grid pattern is given by 

While it is possible to compute the exact power spectrum of such 
a perturbed grid (Gabrielli, 2004, see Appendix C and equation [50] 
below), the general calculation for the Fourier- Taylor expansion is 
very cumbersome. However, three interesting regimes can at least 
partly be discussed qualitatively: 

(i) Very small displacement, = (^^) ^ 1-' in that case, the 
value of s controls the amplitude of the displacement A(i) in equa- 
tion (O and therefore the accuracy of the Fourier- Taylor expansion. 
The smaller s, the better. The worst situation is when s approaches 
diagonal values, e.g., s = (0.5,0.5,0.5). Therefore, to best ex- 
ploit the Fourier- Taylor method on a very slightly perturbed grid, 
it is wise to chose s (and more generally the size of the grid A^g) 
carefully, to minimize as much as possible the amplitude of the dis- 
placements A(i). In the conventions used here, the wisest choice is 
thus s = 0. 

(ii) Significant displacement, ~ 1.- this is the situation of 
Fig. |6] where a ~ 0.88. In that case, there will be always a large 
fraction of particles with large magnitude of the displacement A(i), 
independently of the choice of the offset s. 

(iii) Large displacement, ^ 1: it is difficult in this case to 
make any quantitative statements because the conclusions depend 
on the coherence of the displacement field (see equation |50] be- 
low). However, one can postulate in general that the information 
due to the grid pattern has become subdominant: the particle dis- 
tribution should behave again like the local Poisson realization of 
an isotropic stationary random process and the calculations of § [3] 
should in practice become valid again. 

From this simple discussion, which could be easily extended to the 
more general case A^p 7^ A^|, one sees that the measurement of the 
power spectrum on a perturbed grid has to be performed carefully 
in order to reduce as much as possible the systematic errors on 
the Fourier-Taylor spectrum. However, one expects these errors to 
become significant only when approaching the Nyquist frequency. 

This is well illustrated by Fig. (6) for which we have 
s = (0, 0, 0) and a mean square displacement of order unity [point 
(ii) above]: the magnitude of function g{k) increases quite rapidly 



when k/kny approaches unity. The solid line on each panel still 
represents the function 'y%{k,0). The dotted line corresponds to 
the right member of equation i46\ . but with Pmax = l/A'^p^This 
gives a good idea of the overall behavior of g{k), and this is not 
very surprising. Indeed the calculations of Appendix C (see also, 
e.g., Gabrielli, 2004) give 



[1 - P{l)] 



(50) 



where — 1 ^ p{q) ^ 1 is the congelation coefficient of the dis- 
placement field. From this equation, one sees that for a moderately 
perturbed grid (a^ ^ 1), the power spectrum presents a peak at 
twice the Nyquist frequency of the grid. The amplitude of this peak 
is controlled by the amplitude of the displacements, i.e., the value 
of a. However, when fc^ ^ 1 /a^, the sum in equation l |50| l is dom- 
inated by the g = term, hence 



1 

n: 



k^ > l/a 



(51) 



This means that at wavenumbers large enough, the rough power 
spectrum of the perturbed grid is dominated by the shot noise of 
the particles. Since a ~ 1 in our experiment, this should happen 
for k^ approaching a few units (except for the peak just mentioned 
above at twice the Nyquist frequency of the grid). As a result, the 
effects of the aliasing of the power spectrum should be roughly of 
the order of the shot noise, which is indeed the case on Fig. (6] the 
symbols follow roughly the dotted curve, except perhaps for the 
nearest grid point case (upper left panel). 

Despite the much stronger effects of aliasing on the Fourier- 
Taylor spectrum for the perturbed grid than for the locally Poisso- 
nian case, the systematic errors on P'^'(fc) still decrease rapidly 
with the order, A'^. They remain quite moderate when k/kny is kept 
small enough, for instance k/kny ^1/2 for third order. A'' = 3, 
and increase rapidly if k approaches the Nyquist frequency. 



4.4 Unbiased estimator: a zoom in the range [0, fcny/2] 

From equation i3H we can write an angle-averaged estimator of 



the true power spectrum: 



:,(JV) 



(k) 



P'^'(fc) _ fliv(fc) + l 



(52) 



angles 



Recall that the second term of the right side of this equation has to 
be ignored if we consider a perturbed regular pattern such as initial 
conditions of a A^-body simulation. We explicitly write the angular 
average in this expression, to show how it should be performed to 
get the best estimate of P{k). Equation l l52t gives the estimator we 
propose in that paper. 

In the next section, we shall present a procedure to extend the 
dynamic range of available values of k, which has been restricted 
up to now the the Nyquist range of the (arbitrary) sampling grid 
used to perform the Fourier-Taylor transform. This method will al- 
low us to control the accuracy of the measurements at all wavenum- 
bers, except of course those corresponding to modes close to the 
box size. The previous analyses suggest keeping k far away enough 



^ and of course with P{k) replaced by ProughCi^) (with no shot noise 
correction). 
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5 ARBITRARY WAVENUMBERS WITH CONTROLLED 
ERROR 

In this section, we show how the dynamical range of available val- 
ues of k can be increased arbitrarily while keeping the statistical 
errors and the aliasing effects on the measurements bounded (Jenk- 
ins et al., 1998). 

As shown above, the error on the Fourier-Taylor expansion is 
basically controlled by the magnitude of k: this latter has to remain 
a small enough fraction of the Nyquist frequency to avoid uncon- 
trollable effects of aliasing and hence to be able to maintain suf- 
ficiently low order A'^. Assuming that Sk was computed in a given 
range of values of (each coordinate of) k, [— afcny, afcny], a ^ 1, 
one can now consider an arbitrary shift, ks, such that k + ks is 
outside the available range, e.g., ks — 2akny to have access to the 



range [akny , Sakny 



in ID. Equation ^ is rewritten 



3fc + fes — 



Mi exp(/fc ■ Xi) exp(/fcs 



^^A/;exp(/fc 



Xi) 



(53) 



(54) 
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Figure 7. A zoom in the region [0, fcny/2] of the relative residual, 
5P{k)/P(k) = [Pj:^^ - P{k)]/P{k). on the unbiased third-order power 
spectrum, as measured in our CDM GADGET sample. The upper panel cor- 
responds to present time, where shot-noise correction was performed, while 
the lower one coiresponds to initial conditions, with no shot noise correc- 
tion. 



from kny. The choice proposed in this paper isk ^ fcny/2. Figure|7] 
shows the relative residual on the third-order unbiased estimator in 
the region k G [0, fcny/2]. Such a residual is well below the statis- 
tical errors, of the order of 1.5 x 10"'* and 4 x 10""^ at most in the 
final and the initial stage of the simulation respectively. This error, 
already small, can be reduced further by increasing the order A^. 
However, it does not make sense to do so as long as the statistical 
errors are dominant: to reduce the statistical errors one has to sam- 
ple more modes, i.e., to increase the size of the grid used to perform 
the Fourier-Taylor transform. This is discussed in next section. 



M- = Mi exp(/fcs 



(55) 



One thus obtains, by simple multiplication of the weights by 
exp(/fcs.a;i), a modification of the Fourier- Taylor algorithm that 
gives access to arbitrary values of fc while maintaining the errors 
bounded. Of course, a new set of A'fft transforms has to be per- 
formed for each value of ks . 

While this method allows us to investigate arbitrary values of 
fc, there is still the limitation imposed by the available computer 
resources. For instance increasing the computing volume from a 
cube of sides [— fcny,fcny] to a cube of sides [—M fcny,M fcny] 
amounts to a calculation M^ times more expensive than for the 
original data cube. Instead, one can notice the following property 
of Fourier transform in one dimension. Writing 



52k = ^ ex.p{2Ikxi) 



(56) 



^ exp(2Jfca;i) + 

^ i,Xie[0,L/2l 

+ ^ ejip[2Ik{x,-L/2)+IkL], (57) 



ie[L/2,L[ 



we can set 



n = 2xiiorxi€[0,L/2[, 

Ti = 2xi - Lior Xi £ [L/2,L[. 

We then have 



-^^exp(/fcn), 



(58) 



(59) 



if we assume that kL/{2iT) is an integer; in that case exp{IkL) = 
1. Similarly one can write 



"2fc+l 



1 = ■^Ys{i)exp{Ikri) 



where the sign function S{i) satisfies 

Sii) = 1, if Xi e [0,L/2[, 
S{i) = -1, if Xi e [L/2,L[. 



(60) 



(61) 
(62) 
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Table 3. The sampled values of kL / {2it) when using the folding procedure 
explained in the main text on a grid with Ng = 16. The trivial fundamen- 
tal case, fc = 0, is not shown here. From the left column to the right one, 
one considers M = folding to M = 4 foldings. Each column defines a 
range of available values of fc. As discussed in the text, there can be several 
estimates of P{k) for a given value of k. For instance we have 3 estimates 
of P(4). To make unknown aliasing effects always negligible in practice 
compared to statistical eiTors, we impose k ^ k\^^^ /2, where feiy^' is the 
effective wavenumber probed by the Nyquist frequency of the grid after M 
foldings (this is indicated by the mid-horizontal line on the table). To mini- 
mize the statistical errors, M should then be as small as possible, since this 
parameter controls the sparseness of sampling in Fourier space. To measure 
P(4) for instance, one takes M = 0. In detail, the sampled values of k are 
highlighted in bold on the table, fc = 1, 2, 4, 6, 8, 12, 16, 24, 32, 48, 64; 
Fourier space is increasingly sparse sampled with k. This sparse sampling 
roughly increases linearly in logarithmic scale, but the width of the bin used 
to measure P(k) increases likewise, keeping the number of sampled inde- 
pendent modes C(k) bounded between two values Ci and C2 independent 
of the number of foldings A/, except for M = [here Ci = 98 and 
C2 = 210 for D = 3]. As a result, the statistical errors on the power 
spectrum measurement remain bounded as well, as illustrated by Fig.|9] 



M = 


M = 1 


A/ = 2 


Af = 3 


Af = 4 


C{k) (D = 3) 
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2 
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16 
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2 


4 
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12 


24 


48 


98 


4 


8 


16 


32 


64 


210 


5 


10 


20 


40 


80 


350 


6 


12 


24 


48 


96 


450 


7 


14 


28 


56 


112 


602 


8 


16 


32 


64 


128 


762 



This means that we now have access to a doubly-large range of 
integer values of kL/{2n) by applying Fourier- Taylor method on 
simple foldings of the particle distribution. Note that periodicity is 
not assumed here, except that only the integer values of kL/{2n) 
are available. Such a set is complete in the periodic case. 

Obviously this folding trick can be generalized to higher num- 
ber of dimensions. In that case, the number of Fourier transforms 
needed to perform the calculations increases by a factor 2^ each 
time a factor 2 is gained in the dynamic range of available values 
of fc, exactly as in equations i54\ and J55t . as expected. However, 
we can restrict here to the simplest folding given by equation l |59t 
and ignore foldings involving equation l |60t . This means that we 
are sparse-sampling Fourier space, increasing the dynamic range 
by a factor two each time, maintaining the computational time of 
the same order at each step of the procedure|f]Table[3]shows the list 
of values of fc sampled in one dimension in a simple case, where 
A'^g — 16. Note that the number of available modes per sampled 
value of k remains the same as for the original sampling: the error 
on the rough power spectrum shall remain bounded as we discuss 
later below. 

In practice the folding algorithm works as follows: 

(i) First apply the Fourier- Taylor estimator on the original parti- 
cle distribution sampled on the periodic grid of size A'^g, and mea- 
sure -PeJt ^ (fc) up to some fraction of the Nyquist frequency, afcny. 



° Some of the possible caveats of such a sparse sampling are discussed in 
Jenkins et al. (1998). 



with e.g., a = 1/2 as advocated in the previous section. That corre- 
sponds to the first step of the algorithm, with zero folding, M = 0. 

(ii) (Assuming we are at step M of the process.) Fold the parti- 
cles in each dimension according to equation ( |58t . resample them 
again on a periodic grid of size A'g, but which probes a physical 
size twice smaller than in the previous step, Lm+i = Lm /"i = 
Lo/2^^+^ (Lo = L). Measure P^^^\k) up to a times the Nyquist 
frequency of that grid, which in practice corresponds to twice the 
Nyquist frequency of the grid of the previous step, fc^y^^^^ = 

oi.(M) _ oA-f+ll. 

(iii) Repeat step (ii), until the value a.k[^'^^^ is as large as re- 
quired, for instance until the softening scale of the simulation has 
been reached. 

We now have a set of values of measurements of the power spec- 
trum for a number of ranges of values of fc as for instance illustrated 
by Table |3] These ranges overlap from one folding to another: this 
can be used to check for systematic errors on the measurements. 
Indeed, for one value of fc, there can be several measurements of 
P{k). One has, for each folding, to choose the range of values of fc 
that contribute to the final measurement. To do that, one could think 
of compromising between the statistical errors, which are larger 
when fc is small, and the unknown systematics brought by alias- 
ing, which can become significant when fc approaches the Nyquist 
frequency of the grid. The choice of compromise depends on the or- 
der A'^ of the Fourier- Taylor transform considered. It is theoretically 
possible to find an "optimal" compromise between the parameter a, 
the order A^^ and the resolution of the grid, A'g, to maintain the er- 
rors on the measurement of P(fc) below a given limit at a minimum 
computational cost. However, such a project would go beyond the 
scope of this paper. Our strategy here is rather to make sure that the 
unknown systematics due to aliasing are negligible compared to the 
statistical error given by e.g., equation J43t . We therefore advocate 
a = 1/2 and sufficiently high order Fourier-Taylor approximation. 
For instance, the third-order transform should in practice do well 
enough for A'g of the order or smaller than a thousand. 

Figure [8] illustrates how the measurements behave over the 
range fcL/(27r) G [1, 8192] when one changes the grid resolution 
from A'g = 64 to A'g = 512 in the folding algorithm. In order 
to see the details better, the function represented on each panel is 
^os^t' (fc) /PpD (fc), where Ppd (fc) is the function given by the non- 
linear ansatz of Peacock & Dodds (1996). This figure is supple- 
mented with Fig.|9] which gives for each case considered in Fig. [8] 
the errors on P^'^l (fc), in the following way to make the plot read- 
able: the thin curve (alternatively dotted and continuous) and the 
thick grey curves on the upper part of Fig. |9] correspond to equa- 
tions l|40) and l |43l l respectively. However these equations give the 
error on the estimate of the power spectrum plus the shot noise con- 
tribution of the particles. In order to get an estimate of the statistical 
error on P{k), we compute an expression similar to equation i44t . 
but as follows: 



^ = {E{k) oTEaik)} X 



PpD(fc) + l/A'p 
PpD(fc) 



(63) 



That way, the plot is readable because we replace a noisy func- 
tion with its smooth guess PpD(fc)0 The lower part of the same 



If the error on the rough power spectrum is given by 
2 



AP, 



rough 



rough 



(64) 
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kL/27r kL/277 

Figure 8. The power spectmrn measured at present time in our CDM GADGET sample, witli the estimator Pj^^^ (fe) (eguation l52K using the folding algorithm 
explained in the main text to (sparsely) probe the full dynamic range k = [1, 8192]. To see the details better, the measurements have been divided by the 
analytical proxy PpoC^) of Peacock & Dodds (1996). Each panel corresponds to a choice of the grid resolution, N^, used to perform the Fourier-Taylor 
algorithm in combination with the folding of the particle distribution. The measurements are represented by the thin curve, while the shaded region gives the 
statistical error estimated with equation )40t . The dashed curve corresponds to the shot noise of the particles. The dotted vertical line indicates the softening 
length e: above that scale, P{k) should present a cut-off, which is clearly visible when the signal-to-noise ratio is large enough, see e.g., the lower-right panel. 



figure shows the residual function due to aliasing, Ri,{k){l + 
l/[fpD(fc)Ai'p]} (thick grey, dotted, dashed and sohd hnes, which 
correspond to A'^g = 64, 128, 256 and 512, respectively), as well 



then it follows that the eiTor on the shot-noise-corrected power spectrum 
reads as in equation \AA\ but with the term 1/C(k) replaced with _E'^. To 
obtain a smooth estimate of the errors plotted on Fig.|9] our choice is to re- 
place P{k) in equation )44t with its theoretical proxy, PpY){k). After sim- 
ple algebraic calculations, one just obtains equation (63). In that framework, 
the expressions for the theoretical and measured residuals on the power- 
spectrum estimator can also be written naturally as explained in the main 
text and plotted on Fig.|9] 



as its measurement for A'^g = 128 (symbols). In the last case, 
the function displayed on Fig.|9]is [P^l}{k) - P{k)]/ Ppu{k). It 
should roughly follow the dotted line, which is indeed the case. 
Note that a residual f{k) = J?3(fc){l + l/[PpB{k)Np]} assumes 
= PpD(fc). Since Ppn{k) is a decreasing function of k, the 
function /(fc) is expected to overestimate the true residual^ On 
the contrary, the symbols tend to lie above the "theory", but this 



^'^ We assume here the framework of the assumptions of § |3] For a per- 
turbed grid, as studied in i) |4.3| the residual are expected to be much larger, 
as illustrated by the lower panel of Fig.|7] 
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kL/27T 

Figure 9. The statistical enor on the measurement of P(k) when using the folding procedure explained in the main text. The upper curves con'espond to the 
estimated statistical relative eiTor on P^^^ (fc) (but these curves would not change significantly for other values of the order of the Fourier-Taylor expansion): 
for various values of A^g as indicated on the plot, the thin (alternatively dotted and soHd) and the thick grey curves coiTespond respectively to eqs. )40t and 
)43) for estimating the eiTors on the measurement of P(fe) + 1/Np. The function represented here is given in fact by equation i63\ as discussed in detail in the 
main text. For comparison, the lower curves give an estimate of the expected residual due to aliasing effects on the power spectrum, while the symbol show 
the measurement of these residuals for Ne- = 128. 



disagreement is clearly within statistical errors, which is what re- 
ally matters for the point made here. 

The mid-range values of k in Fig. |9] show that increasing the 
resolution A'g of the grid by a factor of two improves the over- 
all signal-to-noise ratio by a factor 2^^~^^^'^, because C(k) scales 
roughly like k^^^. This is illustrated as well by Fig.jS] where the 
small fluctuations of the measured spectrum at intermediate values 
of k decrease when TVg increases. At small k, however, the errors 
are independent of A'^g as there are fewer and fewer statistically 
independent modes when one approaches the size of the box, what- 
ever the resolution of the grid used to perform the measurements. 
At large k, one is dominated by the shot noise of the particles: when 
P{k) ^ i/Np, as indicated by the dashed line on Figs.[8]and[9l the 
error on P{k) increases dramatically, to become arbitrarily large 
when P{k) is subject to the hard cut-off due to the softening of the 
forces at scales smaller than e, as indicated by the dotted line on 
the figures. Indeed, in the folding method proposed here, one is ac- 



tually able to control the error on the measurement of the quantity 
P{k) + i/Np, and not on the measurement of P{k). 



6 SUMMARY 

In this paper we presented a method to estimate Fourier modes of 
a particle distribution, based on a Taylor expansion of the trigono- 
metric functions. This idea is inspired from the work of Anderson 
& Dahleh (1996). We paid particular attention to the measurement 
of the power spectrum P{k) when the point distribution is the lo- 
cal Poisson realization of a stationary random field, where explicit 
expressions for the ensemble average of a naive rough estimator of 
P{k) were derived. This allowed us to accurately determine the bi- 
ases induced by discreteness and by the Taylor expansion, which 
can be easily corrected for, and to quantify the effects of aliasing, 
which are controlled by the order, A'^, of the expansion. Our calcu- 
lations show that effects of aliasing decrease quickly with A*', as il- 
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lustrated by Table[2l The analytic calculations were confronted suc- 
cessfully with a cosmological A^-body simulation. We also studied 
how the deviations from local Poisson behavior influence the mea- 
surements, such as in initial conditions of the simulation, which 
correspond to a perturbed grid pattern. We proposed an unbiased 
estimator which is nearly free of aliasing, and which still performs 
well for the perturbed grid. The accuracy of this estimator is thus 
entirely controlled by the statistical errors, which arise from the fi- 
nite number of sampled modes. We also showed how the dynamical 
range in Fourier space could be arbitrarily increased while keeping 
the statistical error bounded, by appropriate foldings of the particle 
distribution, as suggested by Jenkins et al. (1998). Note that, while 
the Fourier-Taylor method was applied here to the power spectrum, 
it can be easily generalized to higher order estimators, for instance 
to measure the bispectrum or the trispectrum of the distribution. 

The Fourier- Taylor module as well as the associated power 
spectrum estimator tool we propose is available as an F90 package, 
powmes, at www .projet-horizon.fr or on request from the 
authors. It works with the GADGET file format. 
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APPENDIX A: ENSEMBLE AVERAGES WITH VARIOUS 
ASSUMPTIONS 

There are subtleties that arise when it comes to performing ensem- 
ble averages. While this has been widely discussed in the literature 
(see, e.g., Peebles 1980), we address the issue briefly again here. 
Given a distribution of TVp particles in a (hyper)cubical volume V , 
which is a discrete realization of an underlying random continuous 
field, p{x) (of average unity), one measures the two quantities 



(AI) 



where / and g are some functions. The question we want to address 
here is how to compute the ensemble average of FG over many 
realizations of the underlying distribution, given some constraints. 
There are two relevant cases to consider: 

(i) The "realistic" case: V is a subvolume of a realization of 
much larger volume. In that case ensemble average allows the num- 
ber of particles A*'p to vary, as well as the average density over the 
volume: the quantity 



P{V) 



(A2) 



is allowed to fluctuate around the mean. 

(ii) The N-body simulation standard case: in that case, A'^p is 
fixed, as well as piV) = 1. 

AI The "realistic" case: unconstrained ensemble average 

Following Peebles (1980), we divide V into infinitesimal cells of 
volume 6V, such that they contain zero or one particle. Let n be 
the number of particles they contain, and n — (n) its ensemble 



average. Let p(n, x) be the probability of having n particles in 
an infinitesimal cell at position x. Then, the local random process 
characterizing the realization of the smooth field in a distribution 
of particles gives, assuming that (p) = 1, 



p{l,x) = n5Vp{x), 
p{0,x) = l — n5Vp{x). 



(A3) 
(A4) 



The sum JAIb can be rewritten over the infinitesimal cells, labelled 
as j and at positions Cj , 



So 



which gives in integral notation 

(F) = n / d°xf{x), 
Jv 

and likewise for G. The product, FG, is then 

FG = ^n]/(cj)g(cj) + ^ njnyf{cj)g{c'j). 
3 j^i' 



{FG) = n / f{x)g{x)d''x + 



(A5) 



(A6) 



(A7) 



(A8) 



+ n' / [l + ax,y)]f{x)g{y)d"xd"y. (A9) 
Jv 

In this equation, we have defined the two-point correlation function 

{p{x)p{y)) = l + ^{x,y). (AlO) 

In this paper, we assume stationarity: (,{x, y) — ^{x — y). Equation 
iA9\ is the basis that we used for the calculation of the ensemble av- 
erage of P(^^(A;) in §[13] and fr om which we derive equation l l24b . 



A2 The A'^-body simulation case: constrained ensemble 
average 

Since A'p is now fixed, the method of infinitesimal cells does not 
work anymore, at least not straightforwardly. However, what still 
remains valid is that the probability of having a particle at position 
X is proportional to p{x). More generally, the probability density 
of having a set of particles at positions {xi, ■ ■ ■ ,xmj,), given the 
realization p{x), is given by 

p(xi, ■ ■ • ,a;jvp) = r^p{xi) ■ ■ ■p(2;jvp), (All) 

remembering that p{V) (equation IA2t is now constrained to be 
unity for each realization of the ensemble average. Then 



(F) 



iVp 

E 

8 = 1 



f{Xi)p{xi, ■ ■ ■ ,XNp)d^Xl ■ ■ ■ d^XNp 



= ^ y" fi^)d°x, 



(A12) 



after performing the integrals and then ensemble averaging, using 
(p) = L As a result we converge again to equation l lA7t . since 
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n = Np/V. However, the calculation of (FG) gives, 



{FG) = ^ J^ f{x)g{x)d''x + 



n , iVp(iVp~l) 
V-2 



X / [l + ax,y)]f{x)g{y)d°xd''y, (A13) 



using the definition l lAlOl l. The second term of this expression dif- 
fers from second term of equation l |A9t . by a factor (A'^p — 1)/A'p. 
In addition, the constraint p{V) — 1 for each realization reads, 
after ensemble averaging. 



^{x,y)d°xd^y = 0. 



(A14) 



These differences impose, in particular, that the fundamental mode, 
P'^'(O) is always exactly unity for an A'^-body simulation in our 
choices of units, unlike equation J24t . which was computed using 
the method explained in Appendix Al. 



APPENDIX B: SOME USEFUL ANALYTIC EXPRESSIONS 

The calculation of the number Hi„{k, Af) (equation [23]( can be per- 
formed easily by using the multinomial theorem: 

Ell} 



qi\ X ■ ■ ■ X qoi 

qi-t h9D=" 

X77qi(fcl,A/l) X ■ ■ ■ X rjq^{kD,MD), 



(Bl) 



where 



1/2 



rj„{k,M)= exp[-/ (fc + 27rA/)A] (/ fcA)"dA (B2) 

J -1/2 

is similar to Kn{k,A4) but is computed on a scalar instead of a 
vector: for D = 1, rj„{k,M) = Hn{k,M). Note that rj„{k,M) 
can be computed using the following recursion: 



VoiKM) = (-l)'*'sin(fc/2)/[(A: + 27rA/)/2], 
Vn{k,M) = -izil!!_(fc/2)"7{7"exp(-/fc/2) 



(B3) 



k + 2itM 
r^expil k/2)} + 



n k 



k + 2-RM 
The final result can be expressed as 



Vn-iiKM). (B4) 



2{~l)'"k"nl 
Vn{k,M) = / X 



{k + 27rM)"-t 



1/2 



X <'sin(fc/2)^y^(fc/2 + ^Af)^' 
I i=0 

cos(fc/2) J2 + 



(B5) 



APPENDIX C: THE POWER SPECTRUM OF A 
PERTURBED GRID 

We consider here the case of a three-dimensional grid pattern of 
particles perturbed by a Gaussian random displacement, which is 
curl-free, stationary and isotropic. 



Isotropy and stationarity imply that the joint probability distri- 
bution of displacements Vi — 'P{qi), P2 = 'P{q2) depends only 
on gi2 = |gi — 92 1 and gives 

27 

C{Vi,V2,qi2) 



(27r)3(l -p2)3/2f^6 

Vl + Vl - 2pVl ■ V2 



X exp 



2cr2(l-p2)/3 



where 

2 /^2 



{Vi} = {V^), p(gi2K = (Pi-p2>, 



(CI) 



(C2) 



are the variance of the displacement field and its correlation func- 
tion, respectively. The Fourier modes of the perturbed grid pattern 
are given by equation l |49l l. The constrained ensemble average of 
the power spectrum estimate (keeping A^p fixed; see Appendix A2) 
is 



91 7^92 



5n{0) + 



1 1 

+ 



x{exp{/fc • [P{qi 



91/92 



P(g2)]}). (C3) 
exp[Jfc • (gi - 52)] X 



rf=*Pid^p2/:(Pi,P2,gi2)exp[/fc ■ (Pi -P2)]. (C4) 



Notice that the offset s has disappeared from this expression, as ex- 
pected. After some algebra, one finds equation dSOt , where the di- 
agonal term has been trivially integrated with the off-diagonal one. 
This expression gives the discrete version of the Zel'dovich power 
spectrum (see Schneider & Bartelmann, 1995, for the continuous 
limit). 
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